Direct numerical solutions of the SIR and SEIR models via the Dirichlet series approach

Compartment models are implemented to understand the dynamic of a system. To analyze the models, a numerical tool is required. This manuscript presents an alternative numerical tool for the SIR and SEIR models. The same idea could be applied to other compartment models. The result starts with transforming the SIR model to an equivalent differential equation. The Dirichlet series satisfying the differential equation leads to an alternative numerical method to obtain the model’s solutions. The derived Dirichlet solution not only matches the numerical solution obtained by the fourth-order Runge-Kutta method (RK-4), but it also carries the long-run behavior of the system. The SIR solutions obtained by the RK-4 method, an approximated analytical solution, and the Dirichlet series approximants are graphically compared. The Dirichlet series approximants order 15 and the RK-4 method are almost perfectly matched with the mean square error less than 2 × 10−5. A specific Dirichlet series is considered in the case of the SEIR model. The process to obtain a numerical solution is done in the similar way. The graphical comparisons of the solutions achieved by the Dirichlet series approximants order 20 and the RK-4 method show that both methods produce almost the same solution. The mean square errors of the Dirichlet series approximants order 20 in this case are less than 1.2 × 10−4.


Introduction
The SIR model is a simple system of nonlinear differential equations that has a rich dynamic. This model is an example of compartment models mainly studied in epidemiology [1][2][3]. The SIR model and its adaptive versions have been used broadly to investigate the dynamic of a considered system. The models are applicable in public health to intervene and control a disease's transmission [1,[3][4][5][6][7][8][9][10][11][12]. Once a compartment model appears, a tool is required to analyze the model. Understanding the flow of the model would lead the researchers to an assumption to control the disease spreading. There are several tools available to handle the compartment models depending on the situation we are dealing with. Some researchers have sometimes used a computer model to assess certain situations [13,14]. Others have implemented a numerical method to investigate the models. Since computer models can be utilized only in some circumstances, a numerical method would be an a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 alternative tool. Some recent studies of the compartment models have focused on the fractional-order compartment models by considering the models under the fractional-order α, where α 2 (0, 1). As α approaches 1, the models become the classical compartment models. The fractional-order SIR model has been demonstrated to have biological significance when α � 0.8 [15].
The Dirichlet series is one of the interesting well-known series in the branch of complex analysis since it is a convergent series. The general Dirichlet series is in the form X 1 i¼0 a i e À l i t ; where 0 = λ 0 < λ 1 < . . . and a i 2 C n f0g. The started term a 0 is the approached value of the series as t goes to infinity. The Dirichlet series solution satisfying some specific differential equations has been studied in [16][17][18]. Investigating the Dirichlet series satisfying the compartment models would lead us to a new method to obtain the models' solution which will be mainly discussed in this paper.
This manuscript considers the classical SIR and SEIR models which are equivalent to the fractional-order α when α approaches one. The existence and uniqueness of the models' solutions are confirmed by a well-known theorem appearing in many differential equations books. However, their exact non-parametric solutions, or what we can call the closed-form formulas have been unknown. The results in [19] have illustrated the exact analytical solution of the SIR model in the parametric form. Although the parametric analytical solution has been derived, a numerical method is still required to obtain the approximate solutions S (t), I(t), R(t) at any time t. A well-known numerical method to obtain an approximate solution of any compartment model is the fourth-order Runge-Kutta method (RK-4) which is an iterative method depending on a given step size and it does not provide the long-run behavior of the solution. An approximate closed-form solution on the other hand may give the long-run behavior with no step size needed. This rises up researchers' interest to derive an approximate closed-form solution of the SIR and its related models. Some approximate closed-form solutions of the SIR model have been studied in [20][21][22][23]. Considering the solution in the form of a power series expression is one simple way to obtain an approximate closed-form solution of a nonlinear system of differential equations. However, the power series solution as a function of time t diverges as t approaches infinity. The power series solution of the SIR model has been studied by [24]. A multistage technique repeating the same order of power series approximation with updated initial conditions is a method to get a numerical solution of the SIR model [25]. Using power series approximation or repeatedly using a fixed order of power series approximation can not eventually reach the long-run behavior of the system due to the divergence of the power series. A power series approach together with a defined approximant is presented in [20]. Padé approximant are the other two numerical methods that can lead the approximated solution to the long-run behavior [26,27]. Considering the Dirichlet series solution would be an alternative path to overcome the limitation mentioned above.

SIR model
In the SIR model, the rates that the proportions of susceptible S, infectious I, and recovered R change over time are described by where the initial conditions are assumed to be S(0) = S 0 > 0, I(0) = I 0 > 0, R(0) = 0, and b; g 2 R þ 0 . Another notation used in this model is SðtÞ: Note that the solution functions S, I, R are single variable real-valued functions that S; I; R : R þ 0 ! ½0; 1�. In addition, the total proportion at any time t is assumed to satisfy S(t) + I (t) + R(t) = 1.
Eqs (1)- (3) can be combined to be in the form of a Bernoulli differential equation presented in [19,21]. Solving the Bernoulli differential equation leads to an approximate analytical solution of the model. Here, we present an approach to derive the approximate analytical solution without transforming Eqs (1)-(3) into the Bernoulli differential equation. A result derived from this approach will be used further.
We first rewrite (1) as When we substitute this equation in (2), we get Differentiating both sides of (4) with respect to t and then equating with (5), we have SS 00 À S 02 À ðbS À gÞSS 0 ¼ 0: The solution of this differential equation can be found by defining the function ϕ as Note that our approach of defining ϕ above is different from the approach presented in [19]. After using the notation of ϕ with (6), we have By solving this first-order differential equation and then using (7), we get which is equivalent to The exact solution of (9) is called the analytical solution. By using power series approximation of ln(S/S 0 ), we obtain an approximate solution called a semi-analytical solution. Since we assumed that S 0 + I 0 = 1 and the term ln S S 0 � � can be approximated by S S 0 À 1, it follows that (9) becomes We note here that the assumption I 0 > 0 implies S 0 < 1. It follows that the terms βS 0 − γ and γ − β cannot be zero at the same time. Now, we consider two cases. Case 1: β = γ. Integrating both sides of 10 gives Case 2: β 6 ¼ γ By integrating both sides of 10, we simply get Once we obtain the approximate solution S(t), we can derive the corresponding solutions of I(t) (using (4) and (8)) and R(t) = 1 − S(t) − I(t) as follows: The set of (S(t), I(t), R(t)) for t � 0 presented in (12)- (14) is called a semi-analytical solution.

Dirichlet series solutions of the SIR model
Let u(t) = ln(S(t)/S 0 ). Then, (8) is equivalent to Since this first-order differential equation has no exact non-parametric analytical solution, we consider its Dirichlet series solution. We first set up some notations for ease of computation. Considering 0 = λ 0 < λ 1 < . . . and a i 2 C n f0g, let and let Then (15) is equivalent to F(u, u 0 ) = 0. By Taylor series expansion we have, where A n ¼ u 2;n ðg À bS 0 e u 1;n Þ þ u 0 2;n À ðu 2;n Þ 2 bS 0 e u 1;n 2 þ � � � ¼ a n e À l n t ðg À bS 0 e a 0 À l n Þ þ terms with higher exponents: The following result provides properties of the Dirichlet series solution of (15) which is equivalent to (18).
Proof. It is obvious that a 0 = ln(S 1 /S 0 ) and is a constant and the first term of A 1 is non-constant, a 1 e À l 1 t ðg À bS 0 e a 0 À l 1 Þ. This implies that 0 ¼ ga 0 À bS 0 e a 0 þ b and l 1 ¼ g À bS 0 e a 0 ¼ g À bS 1 . The former is equivalent to Consider (18) with n = 2, observe that the first term (least exponent) of A 2 is a 2 e À l 2 t ðg À bS 0 e a 0 À l 2 Þ 6 ¼ 0 since a 2 6 ¼ 0 and λ 1 < λ 2 . That means that this term would be canceled by a term in Fðu 1;2 ; u 0 1;2 Þ. Note that the exponents of the terms in Fðu 1;2 ; u 0 1;2 Þ are linear combination of λ 1 . This implies that λ 2 = nλ 1 for some integer n � 2. Since Fðu 1;2 ; u 0 1;2 Þ has the term À bS 1 a 2 1 2 e À 2l 1 t ¼ À bS 1 l 2 e À 2l 1 t 6 ¼ 0, it follows that λ 2 = 2λ 1 and then the sum of coefficients of the terms e À l 2 t and e À 2l 1 t must be zero, i.e., a 2 ¼ À bS 1 l 2 l 1 . Now consider (18) with a fixed n > 2 and assume that λ k = kλ 1 for 1 � k � n − 1. The first term of A n is a n e À l n t ðg À bS 0 e a 0 À l n Þ 6 ¼ 0 since a n 6 ¼ 0 and λ 1 < λ n . This means that the exponent −λ n t must appear in Fðu 1;n ; u 0 1;n Þ. Note that Since the exponent −λ n t must be equal to −jλ 1 t for some j � n, it implies that the term with the exponent −λ n t can be extracted from the term À bS 1 e P nÀ 1 i¼1 a i e À il 1 t , see (21). By the series expansion, the term À bS 1 e P nÀ 1 i¼1 a i e À il 1 t can be written as which contains the term of exponent nλ 1 in the form Hence, λ n = nλ 1 and also a n ¼ À bS 1 l n ðnÀ 1Þl 1 . We have done the proof. Note that Eqs (13), (14) and (20) imply that This means that the solution set of the SIR model, . Consider a Dirichlet series approximant which is an approximate form of the solution u(t), say By the previous theorem, to obtain a n we need to get l n which depends on a 1 , . . ., a n−1 and a 1 ¼ a 0 À P n i¼2 a i . It implies that l i for each i = 3, . . ., N may not be unique and the values of a 1 , . . ., a n can be non-unique complex numbers. This provides many approximate forms of the complex Dirichlet series u(t) which are not applicable to computer programming. We need the series u(t) to be a series of real coefficients. The results in Theorem 1 induce us to consider a specific Dirichlet series solution whose exponents satisfy λ n = nλ 1 for n 2 N. The next results are more applicable for achieving the real Dirichlet series u(t).
a i e À l i t and P 1 i¼0 a i ¼ 0, then the n th derivative of f with respect to t at t = 0 satisfies The proof simply follows from the general Leibniz rule and the mathematical induction.
Proof. By the result in Theorem 1, we have b 0 = −a 0 . By plugging the Dirichlet series u(t) into (15) Eq (22) at t = 0 implies that b 1 = βI 0 . By taking the (n − 1) th derivative both sides of (22) with respect to t and then plugging t = 0 together with the result of Lemma 2, we obtain that This implies the result of this theorem.

Dirichlet series approximants of the basic SIR model
Followed by Theorem 1, it is intuitive to consider the case of λ i = iλ 1 . This implies by setting in Theorem 3 that . . .
which can be written in the matrix form as follows: ; where a 0 = ln(S 1 /S 0 ), S 1 can be derived by (20), and b 1 , b 2 , . . . can be obtained by the result of Theorem 3. Now consider an approximate form of the Dirichlet series u N ðtÞ ¼ P N i¼0 a i e À il 1 t satisfying (15). Then by the previous results, we get and a 1 , . . ., a N can be achieved by solving : The existence of the solution a 1 , a 2 , . . ., a N is due to the existence of the inverse of the Vandermonde matrix [28]. Note here that this matrix equation is solvable when N, a 0 , λ 1 , b 1 , . . ., b N are known. After u N (t) is obtained by solving the matrix equation above, an approximation of the solution S(t), say S N (t), is derived by The other approximations of the solutions I(t) and R(t), say I N (t) and R N (t) respectively, are as follows: It is implied by Theorem 3 that c k can be written in terms of b 1 , . . ., b k . The process to get c k presented in Theorem 3 is so applicable that any mathematical programming can perform the task. The terms b k can be derived after the terms c 0 , . . ., c k−2 were known.

Numerical simulations for the SIR model
In this subsection, we present an algorithm to obtain the Dirichlet series approximants S N (t), I N (t), R N (t). According to the previous results, the process of achieving the approximant of the SIR model is as follows: Step 1: Setting c 0 = 1, then c k , k 2 N, can be written in terms of b 1 , b 2 , . . ., b k by using Step 2: Solving 1 − S 1 + (γ/β)ln(S 1 /S 0 ) = 0 for S 1 , setting a 0 = ln(S 1 /S 0 ), λ 1 = r − βS 1 and b 0 = −a 0 .
Step Step 4: Solving the matrix Eq (23) for a 1 , a 2 , . . ., a N . In this step we get u N ðtÞ ¼ P N i¼0 a i e À il 1 t .

A special Dirichlet series solutions of the SEIR model
In this section, we consider another compartment model consisting of the susceptible population S, the exposed population E, the infectious population I, and the recovered population R called the SEIR model whose system of differential equations is as follows: where the initial conditions are assumed to be S(0) = S 0 > 0, E(0) = E 0 > 0, I(0) = I 0 > 0, R(0) = 0, and b; s; g 2 R þ 0 . Additionally, the total proportion at any time t is assumed to satisfy S(t) + E(t) + I(t) + R(t) = 1.
Let us say that we want an approximated solution R(t) of this model to be a specific Dirichlet series in the form:

PLOS ONE
Direct numerical solutions of the SIR and SEIR models via the Dirichlet series approach By plugging (31) in (30) and (29), it is easy to see that the solutions I(t) and E(t) can be written in a similar Dirichlet form of R(t) and the constant terms of the Dirichlet expansion of I(t) and E(t) are zero. This implies that EðtÞ ¼ 0: We also simply obtain By substituting (30) in (27), we get dS ¼ À bSdR: Solving this separable differential equation leads to By (32) and (33), we have Eq (31) with t = 0 implies that Plugging (31) in (30) and setting t = 0 lead to Now setting D n ≔ γ(σE (n) (0) − γI (n) (0)), where E (n) (0) and I (n) (0) are the n th derivative of E(t) and I(t) at t = 0, respectively, for n 2 N [ f0g. Taking the first derivative of (30) with respect to t and then plugging t = 0 lead us to If we keep taking derivatives till the n th derivative of R(t) with respect to t and plugging t = 0, we should get is derived, other solutions can be achieved as follows: The approximated solution (S N (t), E N (t), I N (t), R N (t)) of the SEIR is obtained by choosing a value of λ 1 .

Numerical simulations for the SEIR model
In this subsection, we present a method to achieve an approximated solution of the SEIR model by using the results of the previous subsection. Given the system of differential Eqs (27)-(30) with known parameters β, σ, γ and initial conditions S 0 , E 0 , I 0 , R 0 , the process to derive the approximated solution is the following: Step 1: Derive D n by using (36)-(40).

Error comparison of the SEIR model
The accuracy of our method can be observed by investigating the mean square error.

Conclusions
The presented method needs no step-size interpolation and also shows the long-run behavior of the approximate solutions; however, the exact parametric solution (SIR case) and the solution from the RK-4 method (both SIR and SEIR cases) cannot give the long-run behavior of the system since their processes need step-size interpolation. The exact parametric analytical solution (SIR case) presented in [19] does not give the closed-form solution as a function of time and it still needs a numerical method to get the SIR solution at any time t. The approximate closed-form Dirichlet series solution can be an alternative direct method to get the approximate solution of the SIR and SEIR models.
The theoretical results in this paper provide a numerical algorithm in order to obtain a numerical solution of the SIR and SEIR models. The derived Dirichlet series approximants obviously preserve the positive property of the SIR and SEIR solutions. The comparisons of the solution graphs shown in Figs 1, 2, 4, and 5 assure that the presented method in this manuscript can be used as an alternative method. The SIR solution in the Dirichlet series form obviously gives the long-run behavior (S 1 , 0, 1 − S 1 ) which depends on the given initial conditions and the model's parameters. As a consequence of the approach to the SIR model, a similar process is used with the SEIR model. The assumed characteristic of the solution R(t) as a real Dirichlet series leads us to the approximated solutions of the SEIR model. The solutions

PLOS ONE
Direct numerical solutions of the SIR and SEIR models via the Dirichlet series approach S N (t), E N (t), I N (t), R N (t)) obtained by our method is almost perfectly matched the solution from the RK-4 method. The solutions also obviously preserve the positive property of the SEIR model and carry the long-run behavior (S 1 , E 1 , I 1 , R 1 ) which depends on the initial conditions and parameters. The method presented in this manuscript is not only useful for obtaining the SIR and SEIR models but it also can be used in many kinds of nonlinear dynamical systems.
Supporting information S1 File. Wolfram mathematica codes. There are no raw data used in this paper. All figures in this manuscript are generated by Wolfram Mathematica posted at https://github.com/ Kiattisak-Prathom/Mathematica-Codes-for-PLOS-ONE. (ZIP)